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The two-dimensional nonlinear problem of steady flow past a body submerged beneath an elastic sheet is considered. 
The mathematical model is based on the velocity potential theory with fully nonlinear boundary conditions on the fluid 
boundary and on the elastic sheet, which are coupled throughout the numerical procedure. The integral hodograph 
method is employed to derive the complex velocity potential of the flow which contains the velocity magnitude on 
the interface in explicit form. The coupled problem has been reduced to a system of nonlinear equations with respect 
to the unknown magnitude of the velocity on the interface, which is solved using a collocation method. Case studies 
are undertaken for both subcritical and supercritical flow regimes. Results for interface shape, bending moment and 
pressure distribution are presented for the wide ranges of Froude numbers and depths of submergence. According 
to the dispersion equation, two waves on the interface may exist. The first, longest wave is that caused by gravity, 
and the second, shorter wave is that caused by the elastic sheet. The obtained solution exhibits strongly nonlinear 
interaction of these waves above the submerged body. It is found that near the critical Froude number, there is a range 


of submergences in which the solution does not converge. 


l. INTRODUCTION 


The problem of interaction between a fluid and an elastic 
boundary is a classical problem of fluid mechanics which is of 
interest in offshore and polar engineering, medicine and many 
industrial applications. In the last two decades, this topic has 
received much attention due to the melting of ice in the Arctic 
regions, opening new routes for ships and new regions for re- 
source exploration!*. Most of the studies are devoted to wave 
propagation along the ice sheet, its response on a moving load, 
and effects of heterogeneous properties of the ice sheet such 
as a floe, polynya, cracks, etc.**. 

The works studying the interaction between the flow 
bounded by the ice sheet and the body submerged in the 
fluid started relatively recently. Das and Mandal’ considered 
oblique wave scattering by a circular cylinder submerged be- 
neath a uniform ice sheet of infinite extent and determined the 
transition and reflection coefficients. To solve the problem, 
they employed the multipole expansion method. Sturova® ap- 
plied the method of matched eigenfunction expansions and 
studied the interaction of a submerged cylinder and an inho- 
mogeneous ice sheet including a floe or polynya. Tkacheva? 
considered oscillations of a cylindrical body submerged in a 
fluid beneath an ice sheet and solved the problem through the 
Wiener-Hopf technique. Savin and Savin!? considered the ice 
perturbation by a dipole submerged in water of infinite depth. 
They applied the complex variable technique and the Fourier 
transform to solve the Laplace equation. Shishmarev at al.!! 
studied the strains in an ice cover of a frozen channel which 
are caused by a three-dimensional dipole moving under the 
ice at a constant speed. Li et al.!? considered a circular cylin- 
der submerged below an ice sheet in water of finite depth. The 
solution method is based on the derived Green function which 
satisfies the boundary conditions on the ice/water interface. 
All the works mentioned above regarding submerged bodies 
are based on linear potential flow theory, and the boundary 
value problem is usually formulated in the frequency domain. 


The nonlinear theory of hydroelasticity is currently under 
development. The unknown shape of the ice/fluid interface 
and its higher-order derivatives, which have to satisfy the dy- 
namic boundary condition, are the main challenge to derive 
analytical solutions or develop computational approaches. As 
the dynamic boundary condition gets more complicates, e.g. 
include gravity, surface tension and/or elasticity of the sheet 
covering the fluid, it increases the level of the mathematical 
challenge which has to be addressed. The simplest form of the 
dynamic boundary condition corresponds to free streamline 
flows for which the velocity magnitude on the free streamline 
is assumed to be constant. This class of flows is well devel- 
oped and presented in classical books by Milne-Thomson!?, 
Birkhoff and Zarantonello!*, and Gurevich}. 


For free-surface flows, gravity leads to an additional term 
in the dynamic boundary condition which relates the veloc- 
ity magnitude and the vertical coordinate of the free surface. 
This kind of problem can be reduced to a singular integrodif- 
ferential equation whose form depends on the solution method 
and the choice of the governing functions. Various forms of 
the integro-differential equation were derived by Forbes and 


Schwartz!®, King and Bloor!”, and Faltinsen and Semenov!®. 


For capillary free-surface flows, the dynamic boundary con- 
dition comprises the curvature of the free surface which in- 
volves the first and second derivatives of the free surface. 
Fewer analytical solutions for purely capillary waves are pre- 
sented in literature. Crapper!? developed a closed-form solu- 
tion for a fluid of infinite depth, and Kinnersley”? extended 
his method to a fluid sheet. Crowdy?! developed a method 
based on complex function theory and retrieved Kinnersley’s 
solution in much simpler form and obtained new solutions 
for steady capillary waves on a fluid annulus. However, the 
extension of the method to fluid/structure interaction prob- 
lems with surface tension seems to be nontrivial. Alterna- 
tively, several numerical methods have been developed to 
solve the capillary and capillary-gravity flows. Schwartz and 
Vanden-Broeck” proposed a method based on a boundary- 
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integral formulation and finite difference approximation of the 
derivatives and applied it to the purely capillary and capillary- 
gravity waves. Vanden-Broeck and Miloh”? proposed numeri- 
cal methods based on the Fourier-series expansion and studied 
steep gravity waves. Later, the method was adopted by Blyth 
and Vanden-Broeck”* and Blyth and Părău?’ to compute non- 
linear capillary waves on fluid sheets of finite thickness. Yoon 
and Semenov? considered a cavity flow past a circular cylin- 
der in the presence of surface tension and derived the solution 
using the integral hodograph method. They derived a singu- 
lar integral equation with respect to the velocity magnitude, 
which is solved by the method of successive approximations. 
The method can be applied to solve problems with a more 
complicated form of the dynamic boundary condition which 
comprises higher-order derivatives of the free surface. How- 
ever, the higher-order derivatives of the interface which appear 
in the dynamic boundary condition results to a higher-order 
hypersingular integral equation. A special numerical treat- 
ment is required to solve this type of integral equation. 


The nonlinear theory of hydroelastic waves, for which the 
dynamic boundary condition gets more complicated, has been 
studied intensively in recent decades with emphasis on waves 
generated by a moving load. Most of the studies are focused 
on the analysis and simulation of hydroelastic waves, which 
account for the nonlinearity of the potential flow and elastic 
sheet deformations. ParauPărău and Dias? derived a forced 
nonlinear Schrodinger equation for ice sheet deflection and 
studied the weakly nonlinear effects. Fully nonlinear com- 
putations based on the boundary integral method were pre- 
sented by Guyenne and Parau*. The nonlinear response of 
an infinite ice sheet in the time domain has been studied by 
Bonnefoy et al.8 using a higher-order spectral method. They 
found that at a critical speed at which the linear response is 
infinite, the nonlinear solution remains bounded. Despite the 
progress in the development of numerical methods to solve 
nonlinear problems of hydroelastic waves, their extension to a 
body submerged in fluid beneath an elastic sheet seems to be 
not easy, since the flow potential has to satisfy an additional 
boundary condition on the body surface. 


In the present paper, we study a fully nonlinear problem 
of the hydroelastic waves generated by a body submerged in 
the fluid beneath an elastic sheet. We shall use the model 
of potential flow with the fully nonlinear kinematic and dy- 
namic boundary conditions on the submerged rigid body and 
the elastic sheet which is modelled using the Cosserat theory 
of hyperelastic shells suggested by Plotnikov and Toland??. 
To solve the nonlinear problem, we adopt the solution for a 
cylindrical body moving beneath a free surface??, which is 
obtained using the integral hodograph method. An expres- 
sion for the flow potential which includes the velocity mag- 
nitude on the free surface in an explicit form has been de- 
rived. This gives the possibility to adopt the solution to in 
the present problem, because the velocity magnitude on the 
interface between the fluid and elastic sheet appears in the 
dynamic boundary condition explicitly. The coupling of the 
elastic sheet and fluid problems is based on the condition of 
the same pressure on the interface, one from flow dynamics 
and the second from elastic sheet equilibrium. 


The derivation of the flow potential which contains in ex- 
plicit form the velocity magnitude on the interface and the 
numerical method to solve the coupled fluid/elastic sheet in- 
teraction problem are presented in Section 2. The extended 
numerical results are discussed in Section 3. 


Il. THEORETICAL ANALYSIS 


We consider a two-dimensional steady flow past a cylindri- 
cal body submerged beneath an elastic sheet which is mod- 
elling the ice cover. The characteristic length of the body is 
L, and the thickness of the sheet is B;. A Cartesian coordi- 
nate system XY is defined with the origin at a point inside 
the body and the X-axis along the velocity direction of the 
incoming flow with a constant speed U. The Y—axis points 
vertically upwards. The fluid is assumed to be inviscid and in- 
compressible, and the flow is irrotational. The elastic sheet is 
modelled by the Cosserat theory of hyperelastic shells”°. The 
submerged rigid body is assumed to have an arbitrary shape 
which can be defined by the slope of the body as a function 
of the arc length coordinate S, or Bp = B,(S). The interface 
between the elastic sheet and the liquid is defined by function 
Y(X). The interactions of the submerged body, flow and the 
elastic sheet may generate waves extending to infinity in both 
upstream and downstream directions. On the other hand, the 
flow is uniform at infinity, Y —> œ and —œ < X < œ, and the 
velocity is constant there. In order to provide the same value 
of the flow velocity at infinity in all directions, we introduce 
damping regions P,P, and 7,7) upstream and downstream, 
respectively, where a term providing the wave damping is 
added in the dynamic boundary condition to provide the same 
velocity magnitude at points P) and 7>, or Vp = Vr2 = U. 
Outside the interval P)7>, the flow velocity on the interface 
V(X) =U including infinity. Thus, the fluid surface has a 
limit Y (x) = H which is defined as the submergence of 
the cylinder measured from the origin of the coordinate sys- 
tem. 

We will derive the complex potential of the flow, W = 
W(Z), with Z = X + iY. For the steady flow, the kinematic 
conditions on the body surface and the interface mean that 
the stream function is constant, or 3{W(Z)} = const., as they 
both are streamlines. The dynamic boundary condition on the 
interface is obtained from the Bernoulli equation assuming 
that the hydrodynamic pressure on the interface is the same 
as the pressure conditioning the bending of the elastic sheet. 


y? U? 
p> + Pay + Pice =P > +p8H + Fay (1) 


where U is the speed of the incoming flow, p is the liquid den- 
sity, V = |dW /dZ| is the magnitude of the complex velocity, 
g is the acceleration due to gravity, H is the depth of submer- 
gence, and P, is the atmospheric pressure. The pressure due 
to the bending of the elastic sheet is*+? 


dK 1 
Pee = DB (T+ 5°) +P (2) 
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FIG. 1. (a) Physical plane and (b) the parameter ¢—plane. 


where D' = E /(12(1 — v )) is the flexural rigidity of the elastic 
sheet, K is the curvature of the interface, B; is the thickness of 
the elastic sheet. Equation (2) corresponds to the assumptions: 
the elastic sheet is inextensible and not prestressed”>. 

Two different Froude numbers can be defined based on the 
characteristic length L or the depth of submergence H, respec- 
tively: 


F ole F; aoe (3) 
E > h = . 

vgb Vs 
Using nondimensionalisation based on U, L and p, we have 
v=V/U,x=X/L,y=Y/L,h=H/L, s = S/L, bi = B;/L, 
and W (Z) = ULw(z). Replacing in equation (1) by (2), the 
dynamic boundary condition (1) takes the form 





2(y—h) K 1 
2: 3 
=1 2D ' 4 
: F? (i S ® 
where 
D' dô 
= —, K = —. 
pgLF? ds 


The angle 6 = arcsin(dy/ds) = m + B is the angle between 
the X—axis and the unit tangential vector T which is opposite 
to the velocity direction B. The normal vector n is directed 
from the liquid region outwards, while along the interface the 
spatial coordinate increases in the direction of the vector T 
such that the liquid region is on the left (see Fig. la). 

Equation (4) contains the velocity magnitude along the in- 
terface, the wave elevation y and its derivatives which will be 
related in the following throughout the derived expression for 
the flow potential. 


A. Hodograph method. 


Finding the function w = w(z) directly is a complicated 
problem since the boundary of the flow region is unknown. 
Instead, Joukovskii?! and Michell** proposed to introduce an 
auxiliary parameter plane, or €—plane, which was typically 


chosen as the upper half-plane. Then, they considered two 
functions, which were the complex potential w and the func- 
tion @ = —1In(dw/dz), both as functions of the parameter vari- 
able €. When w(C) and @(¢) are derived, the velocity and the 
flow region can be obtained in the parameter form as follows: 


dw dw dw 


4 
TO, Osaf E ©) 


where the function z = z(zeta) is called as the mapping func- 
tion. 

The flow region beneath the interface and outside the body 
is a doubly connected domain. A canonical region of a doubly 
connected domain is an annulus. By making a cut connecting 
the external and the internal circles of the annulus, the doubly 
connected region becomes simply connected. As shown in 
Fig. la, O-D4 and OD- are the two sides of the cut which 
could have an arbitrary shape but form a right angle with the 
flow boundary at both the body surface (points O_ and O+) 
and the liquid/ice interface (points D_ and D+). 

The simply connected flow region C-D-0OBAO-D4+C+ 
is then transformed into the rectangular domain 
O_D,CD_BAO, in the parameter plane. An upper 
half-plane or unit circle is usually chosen as the parameter 
plane. However, the flow region may have corner points at 
which the mapping z(¢) will be not conformal. Chaplygin 
(see chapter 1(5) in the book!>) pointed out that there is flex- 
ibility in the choice of the region of the parameter variable. 
It should be composed of straight lines and arcs of circles in 
such a way that by means of image transformations of these 
regions it is possible to cover the whole complex plane in a 
simple manner. 

When solving boundary problems, the shape of the auxil- 
iary parameter region is usually chosen with the aim of ob- 
taining the solution of a problem in the simplest form with 
a minimal number of singular points at which the transforma- 
tion of the parameter region onto the complex potential region, 
w, and the region of the function dw/dz is not conformal. In 
the present case of the doubly connected flow region, the addi- 
tional corner points appear at the intersections of two sides of 
the cut and the flow boundary. In order to provide a conformal 
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mapping at these corner points, we have chosen the rectangle 
as the parameter domain, which also has right angles at points 
O_,O,,D_,D4, e.g. the same angles as in the physical plane. 

When the parameter region is chosen as a half-plane or the 
first quadrant, the polynomial functions are usually used to 
construct the mapping function z(¢). Here, for the rectangular 
domain, the polynomial functions will be replaced by Jacobi’s 
theta functions!*, which are quasi-doubly-periodic functions. 
Due the periodicity, they naturally satisfy the same conditions 
on both sides of the cut. Jacobi’s functions have been used to 
solve free surface problems involving doubly connected flow 
regions in the books!*!>-33, 

We may choose the coordinates of the rectangle vertices 
O_04D_Dy4 as (0,0), (7,0), (a,a7/4) and (0,27/4), re- 
spectively, as shown in Fig. la. Here, T is an imaginary num- 
ber. The horizontal length of the rectangle is then equal to 7, 
and its vertical length is equal to 2|t|/4. 

In the flow region, there are two stagnation points marked 
as A, where two streamlines merge into one, and B, where a 
streamline splits into two branches. Positions of these points 
in parameter plane € =a, € = b as well as the position of point 
C, € =c+21/4 which corresponds to infinity in the physical 
plane. The parameters a, b and c should be determined from 
additional conditions at the solution of the problem. 

The interval 0 < € < m on the real axis corresponds to 
the body boundary. The interval c < § < m, in =at/4 
corresponds to part of the interface D_C_, and the interval 
0<€ <c, in = 21/4 corresponds to the other part of the in- 
terface DiC. It should be noticed that points C_ at x 4 —co 
and C}, for x —> +% in the physical plane have been trans- 
formed to the same point C in the parameter region ¢. 


B. Integral hodograph method:derivation of the governing 
functions dw/dz and dw/d¢. 


At this stage we denote the angle of the velocity direction 
along the body as B,(&) and the velocity magnitude on the 
free surface as v(€). With these notations, we have the fol- 
lowing boundary-value problem for the function of complex 
velocity, dw/dz: 








IEE, O<E<m, n=mt/4 © 
dz 
dw Bo(S), O<€ <a, n =0, 
rE =a (Z) =) -B-z a<E<b, n=0, 
z —B,(§)-2n, b<E<a, n=0. 
) 








Wesom = FE nin), O<im<at/4. 8 
In (7) the argument of complex velocity has the jumps equal 
to —7 at stagnation points A (Ẹ = a) and B (€ = b) due to 
the jump of the velocity direction when passing through the 
stagnation point. The two vertical sides of the rectangle in 
the parameter plane correspond to the two sides of the cut in 


4 


the physical plane. The velocities on both sides of the cut 
are the same and therefore the condition of periodicity can be 
applied on the vertical sides the rectangle. The solution of the 
boundary value problem (6)-(8) can be obtained by applying 
the integral formulae derived in?4, 


aie: irai o(e-€) 
7, =e | Th dé In rer) 
_ i fodinv 01(¢ —€ —at/4) T 
tad, a In Fee) Hie] 


It can be easily verified that for 0 < é < m, n = 0 the ar- 
gument, arg[(dw/dz)¢_¢ n-0] = x (%), while for 0 < § < 7, 
in = 17/4, the modulus |dw/dz|¢—¢ in=ar/4 = V(G), i.e. the 
boundary conditions (6) and (7) are satisfied. The boundary 
condition (8) is satisfied due to periodicity of the function 
0(¢). By substituting the boundary conditions (6) and (7) 
into (9) and evaluating the first integral over the step change 
in the function 7(€) at points € =a and € = b, we ob- 
tain the expression for the complex velocity in the rectangle 
0—,04,D-,D4*, 


dw Vi (E =a) (E =b) E * dBp | 








O v (E-— é) 
= Vp p n 
dz Va(E —a)Ba(G—b)  |mJo dE (E-E) 
i [®dinv, H(E— é —21/4) 
n)a dE "(G-E 7/4) 
where Bo is the angle at point O_ which is zero if point O_ is 
the highest point of the body. The constant vp or the velocity 
magnitude at point D+, is determined by satisfying the veloc- 
ity at infinity, É = c+21/4, which is 1, as it has been chosen 
as the reference velocity, or 


dé 











dé io] (10) 


dw 
— =j 11 
dz (11) 


C=c+m1/4 





For steady flows, the stream function y = 3(w) takes con- 
stant values along the body and the interface. According to 
Chaplygin’s special point method!>, to determine the function 
w =w(C), it is sufficient to analyse all special points where 
the mapping is not conformal. These are the stagnation points 
A(€ =a) and B(Ẹ = b) and point C(Ẹ = c + 71/4) corre- 
sponding to infinity in the w—plane. The order of the function 
w —w(C) at these points can be determined by analysing the 
behaviour of the argument of w(¢) in the vicinity of these 
points. Then, the derivation of the complex potential can be 
obtained in the form**. 


dw _ ,Pi($ —a)d4( —a) 1 (E —b) Bal — 5) 
dé 07 (€ —c—mt/4)0?(€ —c+m1/4) 


Dividing (12) by (10), we obtain the derivative of the map- 
ping function as 


dz K 02 (¢ —a)07 (5 —b) 
dë vp 07(€ —c—mt/4)0?(€ —c+a1/4) 





a) 





(13) 








dé +ißol . 


1 (* Bp, O(S—S) 

r| h erae EE 
i [0 dlnv, 0 (E — é — rT/4) 
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whose integration along the intervals 0 < E <candc<&<az 
at n = mt/4 in the €—plane provides the parts D,C and 
D_C_ of the free surface C_C, in ¢—plane, respectively. 
The parameters a,b,c,t and K, and the functions B,(6) and 
d(Inv)/d& have to be determined from physical considera- 
tions and the kinematic boundary condition on the body sur- 
face and the dynamic boundary conditions on the free surface. 


C. System of equations for parameters a,b,c,t and K. 


At infinity, point C_C;(€ =c+7/4), the velocity ap- 
proaches unit (since this velocity is chosen as the reference 
velocity), and its direction is along the X-axis. Therefore, the 
argument of the complex velocity (10) at point c =c+at/4 
should be equal to zero 


EEN =i (14) 
dz} tto 


The scale factor K is determined by the length S, which is the 
perimeter of the body cross-section 


P i = Sp. (15) 
where 

dsp |dz 

de ab lee 





The free surface on the left and right hand sides at infinity 
has the same value of y—coordinate. This is also equivalent 
to that the stream function y = 3(w) is continuous across 
the cut, or 3(wp_) — 3 (wp, ) = 0. By integrating 3(dw/d¢) 
along D_D. passing the point ¢c along a semi-circle C’ of 
an infinitesimal radius €, at which dw/d€ in Eq.(12) has the 
second order singularity, we have 


a(S Taça ATEA Trag) 


z ad cœ d¢ c-e dC 
dw dw 
=3 —dć |) =$ | ix Res — 16 
( œ dọ ) dad ae 


fog gee, } 


Here the first and third terms on the left hand are zero because 
3(w) = const. on the free surface. From this equation it fol- 
lows 


a+b=2ce. (17) 


The depth of submergence, h, and the flowrate, Q, between 
the body and the free surface are related. Therefore, instead of 
a condition for the depth h, we can use the following condition 
for the given flowrate Q, which is the integral of the derivative 
of the complex potential along the side O_D,. 


mt/4 dw 
3 (/ Szat) =Q. (18) 


We may place a vortex with circulation T at the centre of the 
cylinder, which can be nondimensionalized as y = T / (2U L). 
For a circular cylinder, this does not affect the impermeable 
body surface boundary condition, but does change the posi- 
tions of the stagnation points and also affects the free surface 
boundary. For a hydrofoil, y should be determined through 
the Kutta condition at the trailing edge. 

Integrating dw/d¢ along the body surface in the parameter 
plane, we have 


T dw 
x ( 4 Trag) =277. (19) 


In the case y Æ 0, the real part of the potential, @ = R(w), have 
a jump on the sides O_D_ and OD. of the cut, while the 
complex velocity, dw/dz and the stream function y = 3 (w) 
are still continuous across the cut. 

Equations (14) - (19) allow us to determine the unknown 
parameters a,b,c,t and K, which appear in the governing 
equations (10), (12) and (13), once the functions v(&) and 
B,(€) are specified. 


D. Kinematic boundary conditions on the body surface. 


By integrating the modulus of the derivative of the mapping 
function (13) along the side O- O+ in the parameter plane, we 
can obtain the spatial coordinate along the body as a function 
of the parameter variable 


am $ dsp 1 
sy(§) = p gents» (20) 


where dsp/d& = |dz/d¢|c_¢ ņ=0. Since the function p(sp) 
is known, the function B;(&) can be determined from the fol- 
lowing integro-differential equation: 

api _ dBr dsp 


dE 7 dy dE (21) 


By substituting dz/d¢ from (13), this equation takes the form 











dbp K 07 (é —a) 6; (é —b) 
7 Kl) az(E =o 11/4) O2(E —c+nt/4) | 
1 (dpe (E-E) Ve, 
xew] nly ae” BEENS 





i fodinv, 0(§—§'—21/4) |, 
Tı dé! nea (22) 


where K(sp) = dpp /dsp is the curvature of the body. 


E. Nonlinear dynamic boundary condition. 


The dynamic boundary condition (4) includes the interface 
elevation y(s) and its derivatives. Each branch of the interface, 
C_D_(c < § < m) and C,D,(0 < & < c), is evaluated by 
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integration of the derivative of the mapping function (13). The 
parameter form of the interface is as follows, 











5(§ ){cce<n0<E<c} =e wet, (23) 
VE )fccé<m0<é<c} =YDt+S i ue is) 
D {2,0} dË |¢-ë+rr/4 
(24) 
where 
ds _ |dz 
dé 7 at C=€+27/4 
a pee here (25) 
v(5) oF (E —c) 97 lE — c) 


and yp is the vertical coordinate of points D_D., and can be 


obtained from 
z|t|/4 q 
yo =3 {i { Z|) a). 26) 
9 dé | ein 
dB = dB „ds 


The curvature of the interface is 
R 
ds dé’ dé 


where dB/d& is determined by taking the argument of the 
complex velocity from equation (10) at € = E +a71/4, 


7 dw\ 1 f?dinv, | (§ -6&') 
Be) =a8 (Gr) =—5 E rle 





K[s(6)] = (27) 


Pi (5), 
(28) 











2 D (É -a + 77/4) (E -b+ 21/4) 
Pi (E) = —Bo+4 sfin SE 


ap O (E —6'+27/4) | er 
t ah ae fr E G (29) 


and differentiating it respect to variable €, 


dB eee di(E -€') 
de adn db! (aE (EE) a 


Here, prime denotes the derivatives of the functions with re- 
spect to č. The integrand of the above equation has a first- 
order singularity at point 6’ = &, since H (é —€’)~ E —€'. 

The derivatives of the curvature, dk/ds and d?«/ds, can 
be obtained by differentiating (27). They will include higher- 
order derivatives of the function B(&) and a higher-order sin- 
gularity in the integrands, respectively. By substituting the 
derivatives of the curvature into the dynamic boundary con- 
dition (4), we can derive a hypersingular integral equation in 
terms of the function dInv/d6, the solution of which requires 
special treatment. Instead of that, we will determine the func- 
tion v(&) numerically using a collocation method. 














) ag P(E). 


F. Numerical method to determine the function v(é). 


If the tentative function v(&) is given and the system of 
equations (14)-(19) and the integrodifferential equation (22) 
are solved, then the interface z = z(€) depends only on the 
given function v(€). We can chose a fixed set of points Ê, 
k =1,K, distributed on the side D- D, of the parameter region 
corresponding to the interface. Then, the nodes sg = 5(E) and 
the coordinates yg = y(E) in the physical plane are obtained 
using Eqs. (23) and (24). The curvature of the interface and its 
derivatives are determined numerically applying spline inter- 
polation of the nodes yx in the intervals (s,_1,5,). We chose a 
fifth-order spline which provides continuous derivatives along 
the interface up to the fourth derivative, 


y(s) = yk +a gls — Sk—1) +... +anx(s — Sk—1)", 
sk1<S<sk, k=1,...,K, n=5. (31) 
The curvature and its derivatives are obtained by differentiat- 


ing Eq. (31): 


n dk HD _ sg (2 
B=aresiny’, K= y f =) a eee, 
1-y? ds (ye) 





The system of nonlinear equations can be obtained by ap- 
plying the dynamic boundary condition (4) at the points są = 
s(&,), which is written in the form 


Gi(V) = epe(V) — che (V) = 0, 


yasayKy. (32) 
where V = (v1,...,vg)! is the vector of unknown velocities 
vz on the interface; the pressure coefficient on the interface 
due to the flow is 





c) =1- v (33) 
and the pressure coefficient determining the bending of the 
elastic sheet is 


2 
che (V) = 2»| (SS z), +3). (34) 


The system of equations (32) is solved using Newton’s 
method.The Jacobian of the system is evaluated numerically 
using the central difference with Avy, = 1075, k = 1,K. At 
each evaluation of the function G;(V), the system of equations 
(14)-(19) and the integrodifferential equation (22) are solved. 
From 5 to 20 iterations are necessary to get convergence of 
the solution. All solutions, say V*, reported here satisfied the 
condition 


HV Nie 107. (35) 


oe 


which is regarded as giving a sufficiently accurate solution of 
the nonlinear equations. However, the inaccuracy within the 
intervals (E. 4 ; E;) is somewhat lower. It will be discussed in 
detail in section 3.1. 

For supercritical flow regimes, the waving interface may 
extend to infinity. However, due to the finite length of the cal- 
culation region and the condition that the flow is uninform at 
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infinity in all directions (upstream, downstream and at infinite 
depth y — —ce), we need to introduce the damping regions 
P,P, upstream and Ti T) downstream. In these regions, we add 
an artificial term in the boundary condition (4) which may be 
treated as an external applied pressure, 


2_ 2(y—h) d 


v 
a eas (36) 





Cp=l—v 


where the damping coefficient Cy increases from 0 at points 
P, and T; to the values Cg, and Cag at points P) and D, re- 
spectively. The length of the damping regions are chosen to 
be about 2A, where Ag is the wave length of the free sur- 
face progressive waves according to the first approximation 
theory’. 


G. Dispersion equation. 


Differentiating Eq. (4) with respect to the arc length co- 
ordinate along the interface and taking into account that the 
slope of the interface 6 = arcsin(dy/ds) is the angle between 
the unit tangential vector T and the x—axis, we obtain 


dB 3 (2 ) dB 


2dlnv _ 


F’v sinô — F°D 





dst t3 ds ds? 


S 





| . B7) 


The above equation without an elastic sheet (D = 0) and at 
small disturbances of the free surface such that sin ô ~ 6 can 
be written as 
dl 27 
y2 nv ae Zr 


ds A 





(38) 


where the wave length A = Ag = 22F? is known from the 
linear theory of gravity waves*>. In the presence of the sheet 
(D #0), the waves of small amplitude approach a sinusoidal 
curve. Therefore, their slope can be written as 


Ô (s) = max COS (# +6) ; 


where Omax is the amplitude, and @ is the phase of the slope 
relative to point D+ at which s = 0. Then, neglecting the 
square of the curvature, i.e. the second term in brackets (37), 


we obtain 
1 2n\* 
aeo (2) (39) 


According to (38), the coefficient at the angle ô is equal to 
—22/A, where À is the wave length of the interface. There- 
fore, the following equation with respect to wave number 
k =2z/A is obtained 


dlnv 
2 pa 
ds ° 








1 
k= a5 + Dk’. (40) 


This equation may have one, two or no roots, which depends 
on the constant D depending on the thickness of the elastic 
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FIG. 2. Wave number vs. Froude number for different thicknesses of 
the elastic: b; = 0 (solid line), 0.05 (dashed line), 0.1 (dotted line), 
0.2 (dot-dashed line), 0.5 (short dashed line). 


sheet and the Froude number F. The latter case corresponds 
to the subcritical flow regime for which the perturbation on 
the interface decays in both directions. The case of one root 
corresponds to the critical Froude number, Fer, for which the 
waves of the same length Aer = 27 /kcr are extended to infin- 
ity in both directions. Differentiating (40) with respect to k 
and equating the result to zero, after some manipulations we 
obtain 


1 1 

37, ( p8L(1— v?) \ 4 64 EB; 3 

ke = V4 | S= Fe =| =a). 
á va( EB? yo" \ 81 pgL(1 — v?) 


(41) 


For F > Fer, Eq. (40) has two roots, kw < kcr due to gravity 
and kice > ker due to the elastic sheet. We note that the Eq. 
(39) is valid along the whole interface, —oo < x < œ, so both 
waves associated with wave numbers kw and kice may appear 
upstream and downstream of the submerged cylinder. We note 
that the depth of submergence of the body does not appear in 
the dispersion equation, and therefore it does not influence the 
wave number. However, we expect that the depth of submer- 
gence influences the wave amplitude, similar to that observed 


for free-surface flows**. 


The roots of equation (40) for different Froude numbers and 
different thicknesses of the elastic sheet are shown in figure 2. 


Without an elastic sheet, the constant D = 0, and Eq. (40) 
becomes k = 1/F? that corresponds to the hyperbola in fig- 
ure 2. For thickness b; > 0, the critical Froude number ap- 
pears as the minimal Froude to which corresponds only one 
root and the vertical orientation of the slope. The larger rela- 
tive thickness of the elastic sheet results in the larger critical 
Froude number. The dispersion equation predicts only possi- 
ble waves, but the contribution of each wave to the shape of 
the interface have to be determined from the solution of the 
nonlinear problem. 
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Ill. RESULTS 
A. Numerical approach. 


In the discrete form, the solution is sought on a fixed set 
of points ;, j = 1,N, distributed along the side 0-04, 0 < 
é <a, =0, and a fixed set of points Ê, i = 1,M distributed 
along the side D-D4, n = 7/4, of the parameter region. The 
points €; are distributed so as to provide a higher density of 
the points s; = s,(;) near stagnation points A(¢ = a) and 
B(¢ = b). The distribution of the points Ê; is chosen such to 
provide a higher density of the points s; = s(Ê) closer to the 
body and the uniform distribution for |s;| > À. 

The number of nodes on the body and the interface are cho- 
sen in the ranges N = 100+ 200 and M = 1000 + 5000, re- 
spectively, based on the requirement to provide at least 80 
nodes within the shorter waves to get convergence and rea- 
sonable accuracy of the solution. The nodes Ee k=1,K, 
used for interpolation of the interface, yz = y(Ex), sz = s(x), 
are chosen on the set of points Ê; such that i = 4k. Then, 
K = M/4 = 250 = 1000 provides 20 nodes within the shorter 
wave length at which the system of the nonlinear equations 
(32) is solved. 

The integrals appearing in Eq. (10) are evaluated based on 
the linear interpolation of the functions B,(&) and Inv(€) on 
the segments (&;—1,6;) and (E,_1,€;), respectively. They are 
evaluated as follows: 


7 dpp v (E -— é) = Ar; ig; 
= [E m (ZEE) a = afr) +O) 


N. (42) 
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where AB); = By(&;) — Bp(6)-1), Alnv; = Inv(£;) —Inv(&i-1), 
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The integrals (44) - (47) are evaluated using the 8—point 
Legendre-Gauss quadrature formula. The error of the solu- 


tion within the intervals of interpolation (€,_ 1, €,) satisfies the 
relation 





ag, (47) 


*)| < 107° 


ako 
Y IGV 
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(48) 
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which is regarded as giving a sufficiently accurate solution of 
the problem. 

The method of successive approximations is adopted to 
solve the integrodifferential equation (22), which in the dis- 
crete form becomes 


AB B E-B ED 
= F= N 
Agj AG; 

(49) 
where the arc length es the body, s(E), is evaluated us- 
ing (20) with (AB, ;/A&;)“ known at the (k)'" iteration. The 
iteration process converges very fast. After 5 to 10 iterations 
the error is below a prescribed tolerance of 10~°. 

The derivative of the mapping function (13) has a second- 
order singularity at point € = c + mT/4. Therefore, points 
Ê, i = 1,M, along the side D_D. in the parameter region 
are distributed within the two intervals c+ €; < Ê <T, i= 
1,M,, and 0 < Ê <c—&, i = Mı +1,M. These intervals 
correspond to parts D_C_ and D.C of the interface C_C in 
the physical plane. The values £; and €2 are chosen to provide 
the required length of the parts D_C_ and DiC. 





B. Convergence study of the numerical method. 


The formulation of the problem allows us to consider the 
free-surface flow around the submerged circular cylinder if 
we chose zero thickness of the elastic sheet. This case has 
been investigated in** using the method of successive approx- 
imations. The results based on the present collocation method 
and that based on the successive approximations are shown 
in Fig. 3. Without elastic sheet, the submerged body gener- 
ates a progressive wave downstream only. The free surface 
upstream for x < —A tends to be parallel to the x—axis. This 


.M. (43) property was used in the model of Semenov & Wu* to de- 


termine a parameter which affects the velocity magnitude in 
the damping region. Both the present numerical method and 
the computational procedures** predict the same shape of the 
free surface in the region sp) < s < sy}. However, the present 
damping model provides the gradual decay of the wave that is 
necessary to couple the pressure coefficients due to the flow 
and the bending of the elastic sheet in the numerical proce- 
dure. It is seen that the truncation of the computational region 
does not affect the shape of the free surface in the interval 
Spl <S < S71, as seen in Fig. 3. 

In order to investigate the effect of truncation in the pres- 
ence of the sheet, we consider an example of supercritical 
flow for the thickness b; = 0.05 and the depth of submer- 
gence 6.35. From Eq. (40), we obtain the critical Froude 
number, Fẹ, = 1.65. For Froude number F = 2, the flow 
is supercritical, so there are two wave numbers ky = 0.256 
and kice = 0.782 determined from the dispersion equation 
(40). These wave numbers correspond to the wave lengths 
Aw = 24.5R and Aice = 8.03R. The wave length corresponding 
to the linear theory without an elastic sheet is Aj = 2 FR. 
The interface is shown in Fig. 4 for two computational re- 
gions —8 < x/Ao < 9 (solid line) and —6 < x/Ap < 7 (dashed 
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FIG. 3. Verification of the numerical procedure comparing the 
shape of the free surface using different methods: present colloca- 
tion method (solid line); successive approximation?4 (dashed line); 
numerical solution*® (symbols). Froude number F = 2.75, depth of 
submergence h = 7.55. 
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FIG. 4. Flexural-gravity waves generated by the submerged circular 
cylinder. Thickness of elastic sheet bj = 0.05, Froude number F = 2, 
critical Froude number Fer = 1.65, depth of submergence h = 6.35. 


line). In the region —5 < x/Ag < 5, where the damping is ab- 
sent (Cg = 0) for both cases, the interfaces overlap. Outside 
this region, the solid lines and the dashed lines start to di- 
verge. These results show that similar to free-surface flows, 
the truncation of the computational region does not affect 
the part of the interface without damping. The computations 
in Fig. 4 are carried out for the length of the damping re- 
gion, xp, —xp2 = Ao, and the damping region downstream, 
X72 —X71 =2Ap. The values Czy = 2 and Cyr = 10 were used, 
and the numbers of nodes are N = 100 and M = 4000. 











xl, 
(b) 


FIG. 5. Free surface shape (a) and pressure coefficient (b) for Froude 
number F = 1.5, Fer = 1.65 and depth of submergence h = 5.37 
(solid line), 5.60 (dashed line), 6.35 (dotted line) and 11.58 (dash- 
dotted line). 


C. Subcritical flows. 


For Froude numbers F < Fer, Eq. (40) has only complex 
roots, which correspond to decaying perturbations of the inter- 
face caused by the submerged cylinder. In Fig. 5, we show the 
interface profiles for the Froude number F = 1.5 and the rela- 
tive thickness of the elastic sheet b; = 0.05 for different depths 
of submergence. The interface shape is symmetric about the 
y—axis. The shape is different from that observed for the free- 
surface flow without an elastic sheet, for which the free sur- 
face is almost flat upstream and exhibits a wave downstream. 
Thus, the elastic sheet supresses the waves downstream and 
perturbs the flow upstream near the cylinder. As the thickness 
of the elastic sheet tends to zero, the critical Froude number 
F, decreases and become smaller than F. In this case, the 
flow becomes supercritical, which drastically changes the in- 
terface shape. It will be studied in the following subsection. 


It is expected that if cylinder is closer to the elastic sheet, 
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FIG. 6. Free surface shape (a) and Pess coefficient (b) for Froude 
number F = 1.5, Fer = 1.65 and depth of submergence h = 5.37 
(solid line), 5.60 (dashed line), 6.35 (dotted line) and 11.58 (dash- 
dotted line). 


or the depth of submergence is smaller, then interaction be- 
tween the cylinder and the elastic sheet is stronger. This is 
observed in Fig. 5. The deflection of the sheet above the 
cylinder exhibits a trough making the gap between the sheet 
and the cylinder smaller. It is found that there is a minimal, or 
critical depth of submergence, her, below which the numeri- 
cal solution cannot be obtained. For h slightly larger than hcr, 
few iterations to solve the system of nonlinear equations (31) 
are required to get the converged solution, while for h slightly 
smaller than Aer, the elastic sheet starts to oscillate, and the 
iterations do not converge. 


The interface profiles corresponding to for different Froude 
numbers are shown in Fig. 6a. It is seen that the critical depth 
becomes larger as the Froude number approaches the critical 
value . The corresponding distributions of the pressure co- 
efficient (the left axis) and the velocity magnitude (the right 
axis) along the interface are shown in Fig. 6b. For the case of 
free-surface flows, there is also a depth of submergence below 
which a steady free-surface flow does not exist. In that case, 
the velocity magnitude at the crest of the waves tends to zero, 
and the free surface shape forms a corner of . The mechanism 
restricting the existence of the steady flow in the presence of 
the elastic sheet is different. As shown in Fig. 6a, the veloc- 
ity magnitude on the interface is much larger than zero. The 
present results show that there is a condition of consistency in 
interaction between the fluid, elastic sheet and the submerged 
cylinder. 
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FIG. 7. Interface shape at different depths of submergence for Froude 
number F = 3 and thickness of elastic sheet b; = 0.05. 


D. Supercritical flows. 


We begin the computational analysis for relatively high 
Froude number, F = 3, or F / Fer = 1.82, and the ice thickness 
b; = 0.05, for which the dispersion equation (40) has two real 
roots. The corresponding wave numbers are kice = 1.13 and 
ky = 0.1113. The second wave number almost coincides with 
that obtained from the linear theory of gravity waves without 
an elastic sheet, kọ = 1/F? = 0.1111. The ratio of the wave 
lengths generated by the submerged cylinder and the elastic 
sheet is Aice /Aw = 10.11. 

The interfaces near the cylinder are shown in Fig. 7 for 
different depths of submergence. The wave generated by the 
elastic sheet clearly seen upstream at the smaller depths, and 
its amplitude decays as the depth of submergence increases. 
For h = 6.13 the wave upstream almost disappears. The wave 
generated by the cylinder downstream is not completely seen 
since its length exceeds the length of the interface shown in 
Fig. 7. The larger length of the interface is shown in Fig. 8. 
For submergence h = 6.13, the interface coincides with that 
obtained without an elastic sheet which is shown in Fig. 8 by 
the dashed line. As the submergence decreases, the amplitude 
of the wave downstream increases and reaches its maximal 
value at h = 2.95, and then it starts to decrease. This feature 
has been studied in*4. As the cylinder approaches the free 
surface, it affects the free surface at smaller distances from the 
cylinder, and the flow tends to be symmetric about the y—axis, 
similar to that for F — o. 

The elastic sheet weakly influences the interface down- 
stream but generates wave upstream. As the cylinder ap- 
proaches the free surface, the bending of the elastic sheet near 
the cylinder increases, as can be seen from Figs. 7 and 8. This 
causes the increase of the amplitude of the wave upstream. 

The bending moment and the pressure coefficient are shown 
in Figs. 9 and 10 for depths of submergence and , respectively. 
Although the wave due to the elastic sheet is invisible (Fig. 
9c), the pressure coefficient and the bending moment oscillate 
at both directions upstream (left axis in Fig. 9b) and down- 
stream of the cylinder (right axis in Fig. 9b). The frequency 
of oscillations at the upstream is caused by the elastic sheet, 
while the frequency of oscillations at the downstream is due 
to gravity. This qualitatively agrees with results based on the 
linear theory!®!?, 

The results for Froude number F = 2.5, or F/F,, = 1.52, 
and ice thickness b; = 0.05 are shown in Figs. 11-14. The 
wave numbers are as follows: kice = 0.97, ky = 0.1606 and 
ko = 0.1600. The ratio of the wave lengths Ay/Aice = 6.04. 
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FIG. 8. Effect of submergence on the interface shape for the circular 
cylinder at Froude number F = 3. The thickness of elastic sheet 
b; = 0.05 (solid lines) and b; = 0 (dashed lines). 
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FIG. 9. Bending moment (a), pressure coefficient (b) and interface 
shape (c) at Froude number F = 3.0. The thickness of elastic sheet 
bi = 0.05 and depth of submergence h = 6.13. 


The interface shapes for different depths of submergence are 
shown in Fig. 11 and 12. They are similar to those in Figs. 
7 and 8 for F = 3.0. However, the pressure coefficient and 
the bending moment shown in Fig. 13a,b exhibit behaviour 
corresponding to superposition of the gravity wave (longer 
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FIG. 10. The same as in Fig. 9 at the depth of submergence h = 1.94. 
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FIG. 11. Interface shape at different depths of submergence for 
Froude number F = 2.5 and thickness of elastic sheet b; = 0.05. 
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FIG. 12. Effect of submergence on the interface shape for the circular 
cylinder at Froude number F = 2.5. 


wave) and the elastic sheet wave (shorter wave). The ampli- 
tude of the bending moment corresponding to the gravity wave 
is larger than that corresponding to the elastic sheet wave. The 
latter becomes largest at the trough of the gravity wave, and it 
almost disappears at the crest, as seen in Figs. 13b and 14b. 
Such behaviour of the bending moment and the pressure co- 
efficient demonstrates the nonlinear interaction of the elastic 
sheet and the flow, which is still invisible for the interface pro- 
file in Figs. 13c and 14c. 

For Froude number F = 2, or F /F., = 1.21, the wave num- 
bers are as follows: kice = 0.782, ky = 0.256 and kọ = 0.250. 
The ratio of the wave lengths A,,/Aice = 3.05. The interface 
shapes for depths of submergence in the range from 4.6 to 

















FIG. 13. Bending moment (a), pressure coefficient (b) and interface 
shape (c) at Froude number F = 2.5 and the submergence h = 4.4. 
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FIG. 14. The same as in Fig. 13 at the depth of submergence h = 
1.94. 
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FIG. 15. Interface shape at different depths of submergence for 


Froude number F = 2.0. 


11.6 are shown in Fig. 15 and 16. 

At the upstream direction, the wave caused by the elastic 
sheet becomes visible even for the relatively large depth of 
submergence, h = 11.6. Its shape is like a sinusoid with con- 
stant amplitude. The amplitude increases as the depth of sub- 
mergence decreases. At depth of submergence h = 4.61, the 
elastic sheet interacts with the flow in such a way that the wave 
due to gravity extends in the upstream direction, and the inter- 
face exhibits superposition of the both waves. 

At the downstream direction, the interface shape differs 
from a wave of constant amplitude. The shape near the body 
corresponds to the superposition of the gravity and elastic 
sheet waves, and then gradually approaches to the pure gravity 
wave. This also occurs for larger submergence but is less vis- 
ible. The contribution of the elastic sheet wave decays down- 
stream, and the interface approaches the pure gravity wave far 
downstream. The pressure coefficient and the bending mo- 
ment along the interface are shown in Fig. 17a,b for h = 11.65 
and in Fig. 18a,b for h = 4.6. They demonstrate interaction 
between the gravity and the elastic sheet waves. The wave 
due to the elastic sheet is dominating near the cylinder. It 
keeps constant amplitude in the upstream direction and grad- 
ually decays in the downstream direction. For x/Ao > 8, the 
amplitude of the bending moment in Fig. 18a caused by the 
elastic sheet becomes smaller than that caused by the gravity 
wave, and the oscillations become qualitatively similar to that 
in Figs. 13a and 14a. It is seen from Fig. 18a that the distance 
at which the bending moment caused by the elastic sheet de- 
cays is much larger than that for F =2.5. The results for 
Froude number F = 1.7, which is quite close to the critical 
Froude number, F'/F., = 1.03, are shown in Figs. 19-23. The 
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FIG. 16. Effect of submergence on the interface shape for the circular 
cylinder at Froude number F = 2.0. 
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FIG. 17. Bending moment (a), pressure coefficient (b) and interface 
shape (c) at Froude number F = 2.0 and the submergence h = 11.6. 
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FIG. 18. The same as in Fig. 17 at the depth of submergence h = 
4.61. 
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FIG. 19. Interface shape at different depths of submergence for 


Froude number F = 1.7. 
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FIG. 20. Effect of submergence on the interface shape for the circular 
cylinder at Froude number F = 1.7. 
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FIG. 21. Bending moment (a), pressure coefficient (b) and interface 
shape (c) at Froude number F = 1.7 and the submergence h = 11.6. 


wave numbers kice and ky approach each other, and the ratio 
of the wave lengths becomes smaller, Ay, /Ajice. The interfaces 
are shown in Figs. 19 and 20 for the depths of submergence 
from 11.6 to the 6.8. By comparing the interfaces at the depth 
h = 11.6 for the Froude numbers 2.5 and 2 in Figs. 12 and 16, 
respectively, we can find that the amplitude of the interface 
wave increases at the region upstream while the amplitude at 
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FIG. 22. The same as in Fig. 21 at the depth of submergence h = 6.8. 
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FIG. 23. Onset of the solution convergence. 


the region downstream of the cylinder decreases. The latter 
agrees with the free-surface gravity flow past the submerged 
circular cylinder*+. The former indicates the larger effect of 
the elastic sheet at the smaller Froude number. 


The shapes of the interface for different depths of submer- 
gence in Fig. 20 are similar each other, but the amplitudes 
are different. In the upstream direction, x/Ap < 0, the shapes 
are periodic with period of about 2A , that corresponds to the 
superposition of two sinusoidal waves with the gravity and 
elastic sheet waves. 


In the downstream direction, x/Ap > 0, the shape of the 
interface is not exactly periodic because the amplitude of the 
elastic sheet wave decays slowly downstream. By comparing 
the shape of the interfaces, behaviour of the bending moment 
and the pressure coefficient for different Froude numbers, we 
can see that the wave due to the elastic sheet decays slower as 
the Froude number approaches the critical value. The similar 
behaviour of the bending moment and the pressure coefficient 
can be seen in Figs. 21 and 23. 


The interaction between the fluid and elastic sheet and be- 
tween the fluid and the submerged body is not always con- 
sistent for the steady flow. This situation has some analogy 
with the free-surface flow past the submerged cylinder, for 
which there is some combination of the Froude number and 
depth of submergence at which the steady solution does not 
exists. As the depth of submergence increases, the region of 
non-existence of the solution decreases and then disappears 
at large enough submergence. There is a maximal possible 
deflection of the free surface at which the dynamic boundary 
condition can be satisfied. In the presence of the elastic sheet, 
the situation becomes more complicated because the dynamic 
boundary condition includes not only the deflection of the in- 
terface but also its derivatives. The term due to gravity in 
Eq. (31) increases as the Froude number decreases. It may 
cause such combination of the deflection, its derivatives and 
the velocity magnitude on the interface that the pressure dis- 
tribution for the fluid (31) and for the elastic sheet (32) cannot 
be the same, or the dynamic boundary condition (4) cannot 
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be satisfied. For a supercritical flow, the fluid forces domi- 
nate the elastic forces, while for subcritical flows vice versa. 
This changes the flow configuration near the critical Froude 
number. The interface shown in Fig. 19 for F = 1.7 and cor- 
responds to the supercritical flow closest to the onset where 
the converged solution can be obtained. It forms a hill over 
the cylinder, while for the subcritical flow in Fig. 6 for , the 
interface forms a trough. The different limits of the flow con- 
figurations for sub- and supercritical regimes indicate an in- 
consistency of the fluid and elastic forces in some range near 
the critical Froude number. As the submergence increases, the 
deflection of the interface decreases as well as the range of 
Froude numbers at which the steady flow and the elastic sheet 
are not consistent. That can be seen in Fig. 23 where the on- 
set of the region of existence of the steady solution is shown 
in the plane of the parameters Froude number vs. depth of 
submergence. 


IV. CONCLUSIONS. 


A fully nonlinear numerical solution for the problem of 
steady gravity flow past a body submerged beneath an elas- 
tic sheet is presented in the form of the nonlinear analytical 
solution for the fluid part of the problem and the nonlinear 
Cosserat plate model applied to the elastic sheet, which are 
coupled throughout the numerical procedure. The solution of 
the fluid part of the problem is based on the integral hodo- 
graph method employed to construct the complex potential of 
the flow and Jacobi’s elliptic theta functions to deal with the 
doubly connected fluid domain. The curvature and higher- 
order derivatives of the fluid boundary involved in the non- 
linear Cosserat plate model have been evaluated using spline 
interpolation. The coupled problem has been reduced to a sys- 
tem of nonlinear equations with respect to the unknown mag- 
nitude of the velocity on the interface, which are solved using 
a collocation method. 

Steady solutions of the full nonlinear problem were com- 
puted for sub- and supercritical regimes. For subcritical 
regimes, the dispersion equations have no real roots. The 
elastic sheet exhibits a most considerable deflection above the 
cylinder, which rapidly decays away. The deflection forms 
a curve symmetric about the Y-axis with a trough above the 
cylinder. The trough becomes larger as the depth of submer- 
gence decreases. These results qualitatively agree with linear 
solutions!®-!, At the same time, the present nonlinear solu- 
tion revealed a critical depth below which the deflection of 
the interface cannot provide balance between the bending and 
hydrodynamic forces in the steady flow. The critical submer- 
gence increases as the Froude number approaches the critical 
Froude number. 

For supercritical regimes, the dispersion equations have two 
positive real roots which correspond to two wave numbers. 
The smaller wave number is closer to that corresponding to the 
gravity wave behind the cylinder without an elastic sheet, and 
the larger wave number appears in the presence of the elas- 
tic sheet. The dispersion equation does not restrict the flow 
regions in which the waves may occur, i.e. upstream or down- 
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stream of the cylinder. From the linear theories!®!?, it was 


found that the gravity wave occurs downstream of the cylin- 
der, while the wave due to the elastic sheet occurs upstream. 
The present nonlinear solution revealed that the waves may 
occur in both directions, but their amplitudes in each direc- 
tion significantly depend on the perturbation of the interface 
and the ratio F/F.,. 

The calculations are presented for the thickness of the elas- 
tic sheet b; = 0.05 that corresponds to the critical Froude num- 
ber Fe, = 1.65. For larger thickness of the sheet, the critical 
Froude number increases. We expect that the flow configura- 
tions for larger thicknesses of the ice sheet will be similar to 
those for b; = 0.05 at the same ratio F /F,,. 

At high Froude number F'/F,, > 1.8 and depth of submer- 
gence h > 6, the interface shape is almost the same as it is 
without the elastic sheet. However, the effect of the elastic 
sheet can be seen in the behaviour of the bending moment and 
the pressure coefficient at the upstream. They oscillate with 
the wave number corresponding to the elastic sheet. The am- 
plitude of oscillations is larger than that corresponding to the 
gravity wave at the downstream. At smaller submergence, the 
perturbation of the interface increases, and it becomes visi- 
ble at the upstream. At very small depth of submergence, the 
elastic sheet starts to affect the whole interface: at the down- 
stream, the amplitude of the gravity wave slightly decreases; 
and at the upstream, the gravity wave is excited in addition to 
the elastic sheet wave. The interface represents superposition 
of both these waves. 

As the Froude number decreases and approaches the critical 
value F}, the wave caused by the elastic sheet can be observed 
in both directions. Its amplitude is constant in the upstream 
direction. In the downstream direction, the contribution of 
the elastic wave to the resulting shape decays. The closer the 
Froude number F to the critical value F,,, the slower decay is 
observed. 

Similar to subcritical regimes, there is a critical submer- 
gence below which the steady supercritical solution cannot be 
obtained. The closer the Froude number to the critical value 
Fer, the larger the critical depth of submergence. This may 
be caused due to the interaction of the two waves with closer 
wave length that may cause a resonance phenomenon. 

The data that support the findings of this study are available 
from the corresponding author upon reasonable request. 


ly, A. Squire, W. H. Robinson, P. J. Langhorne, and T. G. Haskell, "Vehicles 
and aircraft on floating ice," Nature 333, 159-161 (1988). 

2V, A. Squire, R. J. Hosking, A. D. Kerr, and P. J. Langhorne, "Moving 
Loads on Ice Plates", Kluwer, 1996. 

3A. Korobkin, E. I. Părău, and J.-M. Vanden-Broeck, "The mathematical 
challenges and modelling of hydroelasticity," Phil. Trans. R. Soc. Lond. A 
369, 2803-2812 (2011). 

4P, Guyenne and E. I. Părău, "Computations of fully nonlinear hydroelastic 
solitary waves on deep water," J. Fluid Mech. 713, 307-329 (2012). 

5P, Guyenne and E. I. Părău, "Numerical study of solitary wave attenuation 
in a fragmented ice sheet," Phys. Rev. Fluids, 2, 034002 (2017). 

67. F. Li, G. X. Wu and C. Y. Ji, "Wave radiation and diffraction by a circular 
cylinder submerged below an ice sheet with a crack," J. Fluid Mech. 845, 
682-712 (2018). 

7D. Das and B. N. Mandal, "Oblique wave scattering by a circular cylinder 
submerged beneath an ice-cover," Int. J. Eng. Sci. 44, 166-179 (2006). 

81. V. Sturova, "Radiation of waves by a cylinder submerged in water with 
ice floe or polynya," J. Fluid Mech. 784, 373-395 (2015). 


Nonlinear flexural-gravity waves 


°L. A. Tkacheva, "Oscillations of a cylindrical body submerged in a fluid 
with ice cover," J. Appl. Mech. Tech. Phys. 56, 1084-1095 (2015). 

10A, A. Savin and A. S. Savin, "Ice cover perturbation by a dipole in motion 
within a liquid," Fluid Dyn. 47, 139-146 (2012). 

UK, Shishmarey, T. Khabakhpasheva, A. Korobkin, "Ice response to an un- 
derwater body moving in a frozen channel," Applied Ocean Research 91, 
101877 (2019). 

127, F. Li, G. X. Wu, and Y. Y. Shi, "Interaction of uniform current with 
a circular cylinder submerged below an ice sheet," Appl. Ocean Res. 86, 
310-319 (2019). 

IBL, M. Milne-Thomson, "Theoretical Hydrodynamics," 5th Edition (Dover 
Publications, 1968). 

14G. Birkhoff and E. H. Zarantonello, "Jets, Wakes and Cavities", (Academic 
Press, 1957). 

15M. I. Gurevich, Theory of Jets in Ideal Fluids (Academic Press, 1965). 

16L, K. Forbes and L. W. Schwartz, "Free-surface flow over a semicircular 
obstruction," J. Fluid Mech. 114, 299-314 (1982). 

TAC. King and M. I. G. Bloor, "A semi-inverse method for free-surface flow 
over a submerged body," Q. J. Mech. Appl. Maths 42, 183-202 (1989). 

180, M. Faltinsen and Y.A. Semenov, "The effect of gravity and cavitation on 
a hydrofoil near the free surface," J. Fluid Mech. 597, 371-394 (2008). 

19G, D. Crapper, "An exact solution for progressive capillary waves of arbi- 
trary amplitude," J. Fluid Mech. 2, 532-540 (1957). 

20W. Kinnersley, "Exact large amplitude capillary waves on sheets of fluid," 
J. Fluid Mech., 77, 229-241 (1976). 

21D, G. Crowdy, "Exact solutions for steady capillary waves on a fluid annu- 
lus," J. Nonlinear Sci. 9, 615-640 (1999). 

22L, W. Schwartz and J.-M. Vanden-Broeck, "Numerical solution of the exact 
equations for capillary-gravity waves," J. Flud Mech. 95, 119-139 (1979). 

23J.-M. Vanden-Broeck and T. Miloh, "Computations of steep gravity waves 


15 


by a refinement of Davies-Tulin’s approximation," SIAM J. Appl. Math. 55 
(4), 892-903 (1995). 

24M. G. Blyth and J.-M. Vanden-Broeck, "New solutions for capillary waves 
on fluid sheets," J. Fluid Mech. 507, 255-264 (2004). 

25M. G. Blyth, E. I. Părău, and J.-M. Vanden-Broeck, "Hydroelastic waves 
on fluid sheets," J. Fluid Mech. 689, 541-551 (2011). 

26B-S. Yoon and Y. A. Semenov, "Capillary cavity flow past a circular cylin- 
der," Eur. J. Mech. B Fluids 28, 670-676 (2009). 

?7F, I. Părău and F. Dias, "Nonlinear effects in the response of a floating ice 
plate to a moving load, "J. Fluid Mech. 460, 281-305 (2002). 

28E, Bonnefoy, M. H. Meylan, and P. Ferrant, Nonlinear higher-order spectral 
solution for a two-dimensional moving load on ice, J. Fluid Mech. 621, 
215-242 (2009). 

29P, I, Plotnikov and J. F. Toland, "Modelling nonlinear hydroelastic waves," 
Phil. Trans. R. Soc. Lond. A 369, 2942-2956 (2011). 

30yY A. Semenov and G. X. Wu, "Free-surface gravity flow due to a sub- 
merged body in uniform current," J. Fluid Mech. 883, A60 (2020). 

31N, E. Joukovskii, "Modification of Kirchhoff’s method for determination of 
a fluid motion in two directions at a fixed velocity given on the unknown 
streamline", Math. Sbornik. 15 (1), 121-278 (1890). 

32]. H. Michell, "On the theory of free stream lines", Phil. Trans. R. Soc. 
Lond. A 181, 389-431 (1890). 

33J_ S. Terentiev, A. G. Kirschner, and I. N. Uhlman, "The Hydrodynamics of 
Cavitating Flow." (Backbone Publishing Company, 2011). 

34Y., A. Semenov and G. X. Wu, "Free-surface gravity flow due to a sub- 
merged body in uniform current", J. Fluid Mech. 883, A60 (2020). 

35H, Lamb, Hydrodynamics. Cambridge University Press, (1932). 

36D, Scullen and E. O. Tuck, "Nonlinear free-surface flow computations for 
submerged cylinders", J. Ship Res. 39, 185-193 (1995). 

37F, Smith, A. Korobkin, E. Părău, D. Feltham, and V. Squire, "Modelling of 
sea-ice phenomena," Phil. Trans. R. Soc. Lond. A 376 , 20180157 (2018). 


